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ERROR ESTIMATES FOR APPROXIMATIONS OF DISTRIBUTED ORDER TIME 
FRACTIONAL DIFFUSION WITH NONSMOOTH DATA 

BANGTI JIN, RAYTCHO LAZAROV, DONGWOO SHEEN, AND ZHI ZHOU 


Abstract. In this work, we consider the numerical solution of an initial boundary value problem for 
the distributed order time fractional diffusion equation. The model arises in the mathematical mod¬ 
eling of ultra-slow diffusion processes observed in some physical problems, whose solution decays only 
logarithmically as the time t tends to infinity. We develop a space semidiscrete scheme based on the 
standard Galerkin finite element method, and establish error estimates optimal with respect to data 
regularity in L 2 (Q) and norms for both smooth and nonsmooth initial data. Further, we propose 

two fully discrete schemes, based on the Laplace transform and convolution quadrature generated by the 
backward Euler method, respectively, and provide optimal convergence rates in the L 2 (f7) norm, which 
exhibits exponential convergence and first-order convergence in time, respectively. Extensive numerical 
experiments are provided to verify the error estimates for both smooth and nonsmooth initial data, and 
to examine the asymptotic behavior of the solution. 

Keywords: distributed order, time fractional diffusion, Galerkin finite element method, fully discrete 
scheme, Laplace transform, error estimates 


1. Introduction 


We consider an initial-boundary value problem for the following distributed order time fractional 
diffusion equation for u(x,t): 


D 


M 


(i.i) 


u — Au = f in ff T > t > 0, 
u = 0 on dVt T >t> 0, 
u( 0) = v in 


where SI is a bounded convex polygonal domain in K d (d = 1,2,3) with a boundary <9S2, v is a given 
function on S7, and T > 0 is a fixed value. Here, D[^u denotes the distributed order fractional derivative 
of u in time t (with respect to the weight function p.) defined by 


( 1 . 2 ) 


D 




u(t) = r 

Jo 


D“u(t)^(a) da, 


where D“u, 0 < a < 1, denotes the left-sided Caputo fractional derivative of order a with respect to t 
and it is defined by (see, e.g. [20, p. 91]) 

rt 


(1.3) 


D?u{t) = 


1 


r(l-a) Jo 


[ ( t-s ) 

Jo 


a —u(s) ds, 
ds 


where T(-) denotes Euler’s Gamma function defined by r(x) = / 0 °° t x 1 e t dt, for all x > 0. In this 
paper we consider the case that /i G C[0,1] is a nonzero nonnegative weight function with 0 < /i < 1, 

M0)m(1) > 0. 

In the last three decades, fractional calculus has been extensively studied and successfully employed 
to model anomalous diffusion, in which the mean squared variance grows faster (superdiffusion) or slower 
(subdiffusion) than that in a Gaussian process. The subdiffusion model, which is a diffusion equation 
involving a Caputo fractional derivative D“°u of order au G (0,1) in time: 


(1.4) 


D“°u - A u = f 


a t > t > o, 
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is often employed to model subdiffusion processes in which the mean squared variance grows at a sublincar 
(power type) rate, slower than the linear growth in a Gaussian process for normal diffusion. Formally, the 
subdiffusion model (1.4) can be recovered from the distributed order model (1.1) with a singular weight 
/r(a) = 5(a — cko), where S(a — op) is a Dirac delta function at op- Physically, the subdiffusion process 
can be characterized by a unique diffusion exponent (commonly known as Hurst exponent) showing 
the time dependence of the characteristic displacement [5]. In practice, the physical process may not 
possess a unique Hurst exponent, and the distributed order model (1.1) provides a flexible framework for 
describing a host of continuous and nonstationary signals [5, 6, 41]. Problem (1.1) is frequently applied to 
describe ultraslow diffusion, where the mean squared variance grows only logarithmically with time, e.g., 
Sinai model [40]. The distributed-order fractional model arises often in disordered media, and has been 
successfully used in several applications. For example, Caputo [4] proposed the use of the distributed 
order derivative in generalizing the stress-strain relation in dielectrics, and Atanakovic et al. [1] suggested 
a distributed order wave equation as the constitutive relation for viscoelastic materials to describe stress 
relaxation in a rod. 

In recent years, the theoretical study of problem (1.1) has attracted some attention [34, 21, 28, 33, 27, 
11, 22, 13, 2], Kochubei [21] made some early contributions to the rigorous analysis of the model (1.1), 
by constructing fundamental solutions to the problem and establishing their positivity and subordination 
property. Mainardi et al. [28] studied the existence of a solution, asymptotic behavior, and positivity 
etc. for the case /x(0) > 0 and /x(a)da = c > 0. Meerschaert and Scheffler [34] gave a stochastic 
model for ultraslow diffusion, based on random walks with a random waiting time between jumps whose 
probability tail falls off at a logarithmic rate. Meerschaert et al. [33] provided explicit strong solutions 
and their stochastic analogues. Luchko [27] showed a weak maximum principle for the problem. Li et al. 
[22] established sharp asymptotic behavior of the solution for t —> 0 and t —> oo, in the case of continuous 
density [x with fi{ 1) > 0. Jia et al. [13] studied the well-posedness of a Cauchy problem for an abstract 
distributed-order differential equation using a functional calculus approach. Very recently, Bazhlekova 
[2] analyzed problem (1.1) for fx £ C[0,1], /x > 0 and fi(a) /Oon a set of positive measure. 

The solution to the model (1.1) is rarely available in closed form, which necessitates the development 
of efficient numerical schemes, to enable the successful use of the model (1.1) in practice. Despite the 
extensive studies on the simpler subdiffusion model (1.4) (c/. [23, 30, 46, 7, 45, 36, 16, 17] for an 

incomplete list of works on the numerical approximation of the Caputo fractional derivative D“u(t)), 
there are only very few studies [8, 19, 35] on the distributed order model (1.1). Diethelm and Ford [8] 
developed a numerical scheme for distributed order fractional ODEs. It approximates the distributed 
order derivative D [^u(t) by quadrature, leading to a multi-term time fractional ODE, which can then be 
solved by fractional multi-step methods. Error estimates of the approximation were discussed in [8]. Such 
a technique was also employed to solve nonlinear distributed-order fractional ODEs in [19], but without 
any analysis. Just recently, Morgado and Rebelo [35] developed an implicit finite difference method for 
the model (1.1) with a Lipschitz nonlinear source term in one space dimension. The scheme is based on a 
quadrature approximation of D \^u(t) together with the backward finite difference approximation for the 
Caputo derivative D“u(t), and the second-order finite difference approximation in space. The stability 
of the scheme, and a convergence rate 0(h? + r + (da) 2 ) (with h,r and 5a being the mesh size, time 
step size and step size for quadrature rule, respectively) were established under the assumption that the 
solution u is C 2 in time and C 4 in space and the weight function fx(a) is sufficiently regular. In view of 
the limited smoothing property of the solution operator, cf. Theorem 2.1 below, the regularity required 
by the convergence analysis is restrictive, especially for nonsmooth data. To the best of our knowledge, 
the development of robust numerical schemes for the model (1.1) with nonsmooth data and their rigorous 
analysis have not been carried out, despite its immense practical importance, e.g., in solving inverse 
and/or optimal control problems [18]. 

In this work, we develop a Galerkin finite element method (FEM) for problem (1.1) and establish 
optimal (with respect to data regularity) error estimates for both smooth and nonsmooth initial data v. 
The approximation is based on the finite element space Xh of continuous piecewise linear functions over 
a family of shape regular quasi-uniform partitions {Th\o<h<i of the domain Q into d-simplexes, where h 
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is the maximum diameter of the partition. Then the space semidiscrete Galerkin FEM for problem (1.1) 
is given by: find Uh{t) £ Xh such that 

( L5 ) (D t l]u h,x) + a(u h ,x) = {f,x), Vx £ Xh, T > t > 0, u h (0) = v h , 


where (•,•) denotes the L 2 (fi)-inner product, a(u,w) = (Vu, Vw) for u, w € Hq(£1), and Vh £ X h is an 
approximation of the initial data v. Our default choices for Vh are the L 2 (fl)-projection Vh = Ph,v , for 
v £ L 2 (fl), and the Ritz projection Vh = Rh.v, for Av £ L 2 (fl), where A = —A with a homogeneous 
Dirichlet boundary condition. Further, we develop two fully discrete schemes based on the Laplace 
transform and convolution quadrature generated by the backward Euler method, and provide optimal 
error estimates for both space semidiscrete and fully discrete schemes. 

Our main contributions are as follows. First, in Theorem 2.1, we establish the sharp regularity esti¬ 
mates for the solution to problem (1.1). The proof relies essentially on various refined properties of the 
kernel function w(z) defined in (2.4) in Lemmas 2.1-2.3, which also enable one to apply the established 
techniques for analyzing the semidiscrete and fully discrete schemes. Second, in Theorems 3.1 and 3.2, 
we derive the following error estimates for the space semidiscrete Galerkin scheme (1.5) for t £ (0,T]: 


IK*) - «/.(*) |K ( n) + h\Mu(t) - M*))IK ( n } < {^ 2 * 1 ! 


IMU 2 (fi) 


if v £ L 2 (fi), 
if Av £ L 2 (n). 


For initial data v £ L 2 (fl), the estimates deteriorate as the time t approaches 0, with an extra 
factor in comparison with that for the standard diffusion case [42]. Third, we develop a fully discrete 
scheme based on the Laplace transform. It relies on a contour representation of the semidiscrete solution 
with a hyperbolic contour, and trapezoidal quadrature, cf. Theorem 4.1. Specifically, the fully discrete 
solution UN,h{t) with N + 1 quadrature points satisfies the following error bound for t £ (0, T]: 


IK*) - u N , h (t) IIl 2 ^) < 


c T [e ClN + h 2 |tlog^| IMK (n) 
c (e~ ClN + h 2 ) ||Au||i 2 (n) 


if v £ L 2 (H), 
if Av £ L 2 {n). 


Last, we develop a second fully discrete scheme based on convolution quadrature, generated by the 
backward Euler method, and in Theorem 5.3, establish the first-order convergence of the scheme for both 
smooth and nonsmooth initial data. For example, for nonsmooth initial data v £ L 2 (fl), the fully discrete 
solution approximating the continuous solution u(t n ), t n £ (0, T] (on a uniform grid in time with a 
step size r) satisfies the following bound: 

IK*n) - Uh IU 2 (0) <ct(t + h 2 |log ^| i" 1 |M| L 2 (n) . 

It is worth noting that all error estimates are nearly optimal and expressed in terms of the regularity 
of the initial data directly, and fully verified by extensive numerical experiments. Theoretically, these 
results extend our earlier studies [15, 16, 14] on the subdiffusion model (1.4), which contribute to the 
development and rigorous analysis of robust numerical schemes for the distributed order model (1.1). 

The model (1.1) is closely related to parabolic equations with a positive type memory term, for which 
there are many important studies on numerical schemes based on convolution quadrature [25, 26, 7] and 
Laplace transform [10, 24, 38, 39, 31, 44, 32]. For example, Cuesta et al. [7] developed an abstract 
framework for analyzing convolution quadrature generated by the backward Euler method and second- 
order backward difference. It covers also inhomogeneous and nonlinear problems. Their analysis uses the 
generating function, and the Laplace transform involves w(z) = z a , 0 < a < 1. McLean and Thomee [32] 
studied the Laplace transform method for a fractional order model, whose Laplace transform involves 
w(z) = z a , — 1 < a < 1. These interest works have inspired the current work on the model (1.1). 
However, the existing error analysis does not cover directly the model (1.1), due to the general kernel 
function involved, cf. (2.4). Instead, we shall opt for the general strategy outlined in [26], by deriving 
various refined estimates for the kernel function, especially identifying the suitable condition on the weight 
function ji. These estimates are also essential for analyzing the Laplace transform approach. 

The rest of the paper is organized as follows. In Section 2, we recall the solution theory of the 
mathematical model (1.1) following [21, 22], In Section 3, we develop a space semidiscrete Galerkin 
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scheme, and provide optimal error estimates. Two fully discrete schemes, based on the Laplace transform 
and convolution quadrature, are given in Sections 4 and 5, respectively. Finally, to test and to verify the 
convergence theory, we present in Section 6 extensive numerical experiments. Throughout, the notation 
c, with or without a subscript, denotes a generic constant, which may differ at different occurrences, but 
it is always independent of the mesh size h, the number N of quadrature points, and time step size r. 


2. Solution theory 


In this part, we discuss the solution theory of problem (1.1) using the Laplace transform. The stability 
estimates will play an essential role in developing optimal error estimates. We denote by the Laplace 
transform. First we recall the following well known relation [20, Lemma 2.24, p. 98] 


( 2 . 1 ) 


d?u = z a u — z° 


*(0). 


Next we denote by A the operator — A with a homogeneous Dirichlet boundary condition with a domain 
D{A) = n H 2 (Cl). The H 2 (Cl) regularity of the elliptic problem is essential for the error analysis 

below and it follows from the convexity assumption on the domain Cl. It is well known that the operator 
A generates a bounded analytic semigroup of angle 7r/2, i.e., for any 8 £ (7t/2, n) [12, p. 321, Proposition 
C.4.2] 

1 1 


( 2 . 2 ) 


I \(zI + A) 


-1 I 


< 


" |S(z) 


< 


Vz £ Eg, 


I \ z sin(0)| 

where Eg is a sector with the origin excluded, i.e., 

He = (z 6 C : |arg(z)| < 8}, Eg = Eg \ {0}. 

Now it follows from (1.1) and (2.1) that 

(2.3) zw(z)u(z) + Au{z) = w{z)v, 
where the function w(z) is defined by 

(2.4) w(z) = f z a ~ l p{a) da. 

Jo 

By means of the inverse Laplace transform, the solution u(t) can be represented by 

(2.5) u(t) = S(t)v := —r f e zt H(z)vdz, 

2m Jr e , 5 

where the kernel H(z) is defined by 

H(z) = (zw(z)I + Ay'wiz), 

and the contour by 

(2.6) Tg^ = {z £ C : \z\ = 6, \ &rgz\ < 8} U {z £ C : z = pe ±l9 , p > <5}. 

We begin by discussing the regularity estimates of the solution. To this end, we first give a few 
elementary properties of the function w(z). The first is the sector-preserving property, which enables 
applying the resolvent estimate (2.2) in the error analysis to be developed in Sections 3-5. 

Lemma 2.1. Let 8 £ (7t/2,7t) and assume that p(0)p(l) > 0. Then zw(z) £ Eg/ with 8' £ (7t/ 2, 7r) for 
all z £ Eg and 8' depends only on p and 8. 

Proof. Let z = re lv with ip £ \—9,8], If ip £ (—7t/2, 7t/ 2), we have 

3 %(zw(z)) = / r“ cos (aip)p(a) da > 0. 

Jo 
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It suffices to consider the case ip £ [t/2, 9}. First we claim that there exists an ro £ (0,1) only dependent 
on p such that 5R {zw{z)) > 0 for all r < rg. By the assumption p(0) > 0, we can find a small eo > 0 such 
that min Q , g [ 0jeo ] cos(aTr)p(a) = <5q > 0. Hence 


(2.7) 


r£ o /*! 

5ft(zu>(z)) > / r a cos(ap)p(a) da — / r“| cos(ay>)|/x(a:) da 

Jo j Cn 


> 5o / r“ da - ||/Lt||c[o,i] / r“ da 
J 0 J €q 


> -(Inr) 1 [<5 0 — r e °(<5o + ||mIIc[o,i])] • 


Then direct calculation yields 


(5 0 - r e °( 6 0 + |Hlc[o,i]) >0 Vr < r 0 -■ ^ 




l/«o 


e (0,1), 


+ HmI|c[0,i] / 

and the desired claim 5ft(2u;(z)) > 0 follows. Now we consider the case r > ro and ip £ [7t/ 2,0]. In fact, 


|tan(arg(zw;( 2 )))| = — 75 - 


/„ r" sin(a<p)f/,(a) da\ r a sm(anp)p(a) da 


| f n r“ cos(ap)p(a) da | ||Hlc[o,i] f 0 r a da 

In view of the assumption p{\) > 0 and p £ C[0,1], we may find a small ei > 0 such that min Q , g n_ eil ] p(a) > 
<5i > 0 and 

/ r a sm(ap)p(a) da > / r a sin(ap)p(a) da 
JO J l-Ei 

> f r“ sin(0)^i(a) da > 5i sin(0) f r" da. 

J 1-ei J 1-ei 


For r > ro, clearly there holds 

/ r a da = f r; 


( 2 . 8 ) 


' 1—ei 


> T0£l 


11— ei 0 \n> 

/'(- 

/o V r 0 


da > r q 


ll-ei \ r o 

r -1 


da 


da > ro£i / r“ da. 

Jo 


□ 


Then we have 

(2.9) | tan(arg(zw(z)))| > <5i sin(0)r o ei/||/Lt|| C [o,i] = : c'. 

Hence zw{z) £ Eg, with 9' = t: — arctan(c'). 

Remark 2.1. FFe note that the constants 61 , r 0 and e\ in (2.9) are all independent of the choice 9. Hence 
in case of 9 = ir — e for a small e > 0, zw{z) £ I 2$> with 9' = 7r — arctan(csin(0)) = 7r — arctan(csin(e)) 

7t - ce, where c = <Si7- 0 ei/||/i||c[o,i]- 

The second result is an upper bound on the kernel w(z), which can be obtained by an elementary 
calculation. 

Lemma 2.2. Let p £ C[0,1] be a nonnegative function. Then there holds 

The third result gives a lower bound on the function zw(z). 

Lemma 2.3. Let 9 £ ( 71 / 2 , 71 ) and assume that p(0)p(l) > 0. Then there exists a constant c > 0 
dependent only on 9 and p such that for any z £ Eg 

(2.10) 1 ' 11 ' ' ~ ’ |z| “ 1 

(2.11) 


\zw(z) I > c / r a da = c- -—, 

Jo log \z\ 

|z|iu(| 2 |) > \zw(z)\ > c|.z|tu(|,z|). 
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Proof. Let z = re llp . Using p(l) > 0 and /i £ C[0, 1], we can find a small ei > 0 such that min ae r 1 _ £lil i p{a) > 
<5i > 0. Then we have for all r > 1 

[ r a p(a ) da > f r a p{a) da > Si f r a da > ei<5j f r“ da. 

Jo J 1-ex J 1-ei JO 

Similarly, we may find a small e 2 > 0 such that min Q , g [ 0 e2 ] p{a) > 82 > 0 and then for all r < 1 

f r a p(a ) da > 6282 f r a da. 

Jo Jo 

Hence for ip £ (6 — n, tt — 9), we get for Ci = min(ei5i, € 282 ) 

\zw(z)\ > (zw(z))> cos (71 — 9) I r a p(a) da > c\ cos(7r — 0) r a da. 

Jo Jo 

Now it suffices to consider the case <p £ [7 t — 9,9], and the case ip £ [— 9, 9 — tt] follows analogously. From 
(2.7), we deduce 

"1 / J Z \ l/ e O 


\zw(z)\ > $i(zw(z)) > 


r" da Vr < ro = 


2((5 0 + ||/j||c[o,i])7 


Then a similar argument for deriving (2.8) shows that the inequality (2.10) holds for r > rg and <p £ 
[tt — 9,9], thereby showing (2.10). The inequality (2.11) follows from 

IHc[o,i] / T a da > f r a p{a) da = \z\w{\z\) 

Jo Jo 

and the trivial inequality |zta(z)| < |,;|mj(| 2 |). □ 

Now we give the main result of this section, namely, stability of problem (1.1) with / = 0. 

Theorem 2.1. Let p £ C[0,1] be a non-negative function with p(0)p(l) > 0. Then the solution u to 
problem (1.1) with / = 0 satisfies the following stability estimates for t £ (0,T] and v = 0,1: 

(2.12) ||H^( m )(t)u|| L2(n) < c T t- m -'l 1 (ty\\v\\ L , ( n ) , v £ L 2 (H),m > 0, 

(2.13) \\A v S^\t)v\\ L 2 m < ct- m+ 1 - l '£ 2 (t) 1 - v \\Av\\ L 2 {n) , v £ D(A), v + m> 1, 

where t\(t) = (log(2T/t)) _1 , ^(t) = log (max(t^ 1 ,2)) and ct > 0 is a constant that may depend on d, 
H, p, M, to and T. 


Proof. The existence and uniqueness of a weak solution was already shown in [22], and it suffices to show 
the stability estimates (2.12) and (2.13). First, by the resolvent estimate (2.2), we obtain the following 
basic estimate on the kernel H(z ) 

||lL(z)|| = | \(zw(z)I + H) -1 |||wj(z)| < M/\z\ \/z £ Eg. 

Let t > 0, 9 £ ( 7 t/ 2 , 7r), <5 > 0. We choose <5 = 1/t and denote for short T = Tg^. First we derive the 
estimate (2.12) for v = 0 and m > 0. By the solution representation (2.5), we deduce 


\\s (m) m = 


< c 


~~7 / z m e zt H(z) dz < c [ \z\ m e^ z)t \\H(z)\\ \dz\ 

27T1 Jy Jy 

dr + J e cos ^t~ m d'l/j'j <ct~ m . 


jm— 1 rt cos 0 

r e 


Next we prove estimate (2.12) for v = 1 and to > 0. To this end, we take 8 = 2T/f in the contour F. By 
applying the operator A to both sides of (2.5) and differentiating with respect to time t we arrive at 


(2.14) 

Owing to the identity 


AS^it) = -i-r / z m e zt AH(z)dz. 
27ri Jy 


AH(z) = A(zw(z)I + A)~ 1 w{z) = (/ - zH(z))w(z), 
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it follows from Lemma 2.2 that 
(2.15) 

Hence we obtain from (2.14) 
||^ m )(t)||<c/ \z 


||Aff(;z)|| < c|iy(z)| < c-j— 77 — 7-7 Vz e Efl. 


kllog^l 


r 

-OO 


< c 


\z\~l 
kl log |^| 
1-1 - 1 


0 5R(z)i 


l<M 


J 2 T/t iogT 

Since 2 T/t > 2, we can bound the first term I by 

I < 


r- ‘7- e rt cos 9 dr + c r ^~ m , 2r /* , * f e 2Tcos ^ di/j =: I + II. 

\og(2T/t) J_q 


nOO /*00 

/ r m (logr)- 1 e rtcos0 dr < c^i(t) / r "V tcos0 dr < c T t“ m -i 4 (t). 

d2T/t J2T/t 

Meanwhile the second term // can be bounded by 

II = c T t- m {2T/t - l)/log(2T/t) < c T t- m {2T/t)/\og{2T/t) = C 7 d“ m_ 1 £i(t). 


This shows the first estimate (2.12). To prove the second estimate (2.13) with v = 0, we choose 6 = 1/t 
and denote again T = Tg^. Then 

S^ m \t)v= -L [ z m e zt H{z)vdz= -L [ z m ~ 1 e zt zA~ 1 H(z)Avdz. 

2tt 1 J r 27 ti Jp 

Upon noting the identity 

zH _ 1 7d(;z) = zw(z)A~ 1 (zw(z)I + H ) _1 = H -1 - ( zw{z)I + H ) _1 
and the fact that f r z m ~ 1 e zt dz = 0 for m > 1 , we have 

S( m \t)v = ^-- [ 7 / 2 m " 1 e zt (xu)(x)/ + H)- 1 dzHu 

27T1 Jp 27T1 Jp 

,m- 1 e zt(, u; (,) / + A yl dzAv 

By (2.2) and Lemma 2.3 we obtain 

||( 2 «;( 2 )/ + H) _i || < M\zw{z)\~ l < c l0g ^| , 

\z\ - 1 

and thus using this estimate and the the monotone decreasing property of the function f(x) = — 
on the positive real axis R + , we get 



||S (m) (*WU’ ( n) < c Qf \z\ m ~ 1 e^ z ' >t \\(zw(z)I + |dz|) ||Hu|| i2(n) 

<c( e rtcos 9 dr+ [ e cos ^ dip) || Ai 

\Ji/t r-1 1/t — 1 J_g J 


\L 2 (Q) 


< ct ~m+l]^Lll\\ Av \\ L2{ny 


We observe that if t n l > 2, i.e. t n < 1/2, then lu ^" 1 < 21og(t rl 1 ). Otherwise if t n l < 2, i.e. t n > 1/2, 

then by the monotonicity of the function f{x) = on R + , lo ^” 1 = < 21og(2). Then we 

deduce 

||S' (m) (t)u|| i 2 (f 2 ) < ct- m+ H 2 {t)\\Av\\ L ,^). 

Lastly, note that (2.13) with v = 1 is equivalent to (2.12) with v = 0 and v replaced by Av. This 
completes the proof of the theorem. □ 
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Remark 2.2. The a priori estimate of the solution at short time is given in Theorem 2.1, in which the 
constant ct depends on the final time T (see also [22, Theorem 2.2] for the special case v G L 2 (Cl), v = 1 
and m = 0). The long time asymptotic behavior of the solution in case of v G L)(A) was given in [22, 
Theorem 2.1], i.e., it decays like (logt ) -1 as t —>• oo; see also [43, example 6.5] for related discussions on 
asymptotic decay. 


3. Semidiscrete discretization by Galerkin FEM 


Now we discuss the space semidiscrete scheme (1.5) based on the Galerkin FEM. On the finite element 
space Xh, we define the T 2 (fl)-orthogonal projection Ph : L 2 (Li) — > Xh and the Ritz projection Rh : 
H(j (fi) -»• X h , respectively, by 

(Ph<P, X) = (<A X) Vx £ I/,, 

{VR h y, Vx) = (V<p, Vx) Vx g x h . 

The Ritz projection Rh and the L 2 (fl)-projection Ph have the following properties [42]. 

Lemma 3.1. Let the mesh Th be quasi-uniform. Then the operators Rh and Ph satisfy: 

\\<p - Rh<fi\\L*((i) + h\\V(y - RhP)\\mu) < ch q \\ip\\ Hq (ty My G nHi(n), <7=1,2, 

\\<p - ^lU’(n) + h \\v(y - PhT)\\mn) < cW\\y\\ m{n) My g ^(n) n h*(si), q = 1 , 2 . 

In addition, Ph is stable on for 0 < q < 1. 


The space semidiscrete Galerkin scheme for problem (1.1) reads: find Uh{t) G Xh such that 

(3.1) (D Y ] u h , x) + (Vu h , Vx) = (/, x) Vx G X h , 

with 117 ,( 0 ) = Vh G Xh . Upon introducing the discrete Laplacian A h '■ Xh —> Xh defined by 

~{X h y, X) = (V<A Vx) My, x G X h , 

and Ah = —Xh, the space semidiscrete Galerkin scheme (3.1) can be rewritten as 

(3.2) D Y ] u h {t) + A h u h {t) = 0, t > 0 
with u h ( 0) = v h G X h and A h = -A h . 

For the error analysis of the semidiscrete scheme (3.2), we employ an operator trick due to Fujita and 
Suzuki [9]. To this end, we first represent the semidiscrete solution Uh to (3.2) by 

(3.3) u h {t) = S h {t)v h ■= J~r [ e zt (zw(z)I + A h )~ 1 w(z)v h dz. 

2m Jve, s 


Lemma 3.2. For any if G Hq(CI) and z G E for 9 G (7 t/2 , tt), there holds 
(3-4) \zw(z)\\\if\\ 2 L 2 {n) + HVV’Hl^n) < c zw(z)\\ip\\l 2{n) + ||V-0 || 2 


Proof. With Lemma 2.1, the proof is identical to that of [3, Lemma 3.3], and hence omitted. □ 

Now we introduce the error function e(f) := u(t) — Uh(t) which, in view of (2.5) and (3.3), can be 
represented by 

(3.5) e(t) =-[ e zt w(z){fy(z)-y h {z))dz, 

2m he,s 

with y(z) = (zw(z)I + A)~^v and y>h{z) = ( zw(z)I + Ah)~ 1 PhV ■ The following lemma shows a bound 
on the error yh — y. It follows directly from Lemma 3.2, similar to [3, Lemma 3.4], and hence the proof 
is omitted. 


Lemma 3.3. Let v G L 2 (Ll), z G Eg with 9 G ( 7 t/ 2 , 7 t), y(z) = ( zw(z)I + A) x v and yh(z) = (zw(z)I + 
Ah)~ 1 PhV- Then there holds 

(3.6) || y{z) - y h {z)\\ L 2(n) + h\\X 7 (y(z) - y h {z))\\ L 2(n) < ch 2 \\v\\ L 2^ n) . 




ERROR ESTIMATES FOR DISTRIBUTED ORDER FPDE 


9 


Now we can state an error estimate for nonsmooth initial data v G L 2 (H). 

Theorem 3.1. Let u and Uh be the solutions of problem (1.1) and (3.2) with v G L 2 (ft) and Vh = PhV, 
respectively. Then for t > 0 and i\(t) = log(2T/f) _1 , there holds 

IK*) - u h (t) || L 2 (n) + h\\X7(u(t) - u h (t))\\ L 2 {Q) < c T h 2 t" 1 £ 1 (t)||w|| i 2 (n) . 

Proof. In the error representation (3.5), by choosing 6 = 2 T/t in the contour Tqj and appealing to 
Lemmas 3.3 and 2.2, we deduce 

l|Ve(i)||„, 0 , < ch Qs‘™'I^±dr\H L , (n)+ ch :=/ + //. 

Now the first term I can be bounded by 

roo | 7 pOO 

I <ch e rtcose - - dr |M| L 2 (n) < / e rtc ° s8 dr < c T ht-Hfft)\\v\\ L ^ n) , 

J 2 T/t logr log(2 T/t) J 2T/t 

and the second term II is bounded by 

H ~ t\og(2T/t) < CT/it _1 4(t)||w||z,2 (n) . 

The bound on || Ve(£) ||x, 2 (f 2 ) now follows by the triangle inequality. A similar argument yields the desired 
L 2 (Q) error estimate. □ 

Next we turn to the case of smooth initial data, i.e., Av G L 2 (fl ), and derive the following error 
estimate. 


Theorem 3.2. Let u and Uh be the solutions of problem (1.1) and (3.2) with v G H 2 (Ll) and Vh = RhV, 
respectively. Then for t > 0, there holds: 

(3-7) \\u(t) - u h (t) || L 2 (n) + h\\S7(u(t) - u h (t)) || L 2 (n) < ch 2 \\Av\\ L 2 (n) . 

Proof. Like before, we take 9 G (7t/2,7t) and 5 = 1/t in the contour 1^/. Then the error ehit) = 
u(t) — Uh(t) can be represented by 


Using the identity 

we deduce 
(3.8) 


e h{t) = v—7 f e zt w(z) ((zw(z)I + A) 1 - (zw(z)I + A h ) 1 R h )vdz. 

2m jr e .s 

w(z)(zw(z)I + A) -1 = z~ x I - z~ 1 {zw{z)I + A)~ l A, 

= 27 ri iyj e Zt *- 1 (fa( z )-?(z))dz + j e^z^fv - R h v)dz^j 


where (p{z) = (zw(z)I + A) l Av and fih(z) = ( zw(z)I + Ah) l AhRhV. Then Lemmas 3.1 and 3.3, and 
the identity AhRh = Ph A give 

\\v(z) ~ Vh(z)\\ L 2 (n) + h\\V(<p(z) - yh{z))\\ L i(a.) < ch 2 \\Av || L 2 (n) . 

Now it follows from this and the representation (3.8) that 


||e fc (t)|| < c/i 2 ||An|| L 2 (n) ^Je rtcose r- 1 dr + e cos ^ dif j < ch 2 \\Av 
which gives the L 2 (U)-error estimate. The H l (Ll) estimate follows analogously. 




□ 


Remark 3.1. The error estimate for nonsmooth initial data v G L 2 (fl) deteriorates like t~ 1 £i(t) as 
t —> 0 + . The behavior agrees with the solution singularity in Theorem 2.1. The factor t _1 £i(f) is 
different from that for subdiffusion [15] and multi-term time fractional diffusion [14]. In contrast, for 
smooth initial data Av G L 2 (fl), the error estimate is uniform in t. 
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4. Fully discrete scheme I: Laplace transform 


The first fully discrete scheme is based on the Laplace transform. To this end, we select a proper contour 
Tg^ in the integral representation (3.3) of the semidiscrete solution uu, and then apply a quadrature rule. 
We follow the works [ 0, 24, 32, 38, 39, 44] and deform the contour T^ to be a curve with the following 
parametric representation 


(4.1) 


z(0 := A(1 + sin(i£ - ip)), 


with A > 0, ip £ (0, 7 t/ 2 ) and (el. The optimal choices of A and ip will be given below, in the proof 
of Lemma 4.4. This deformation is valid since it does not transverse the poles of the kernel function 
H{z)v = (zw(z) + Ah)^ 1 w(z)v, cf., Lemma 2.1 and Lemma 4.3 below. Upon letting z = x + iy, we 
deduce that the contour (4.1) is the left branch of the hyperbola 


— A 


y 


= L 


(4.2) . , . 

\ A sin ip j \\ cos ip , 

which intersects the real axis at x = A(1 ± sin^i) and has asymptotes y = ±(A — x) cot ip. Now we can 
represent the semidiscrete solution Uh(t) by 


(4.3) 


/ OO 

g(€,t)d£ 

-co 


with the integrand g(t;, t) being defined by 


(4.4) 


<?(£, t) = —e z{i)t (z{£)w(z(£))I + A h ) 1 w(z(£))z'(£)v h . 


Remark 4.1. The integrand <?(£,£) exhibits a double exponential decay as |£| -A oo for t > 0. 


Now we describe the quadrature rule for approximating (4.3). By setting Zj = z(fj) and z' :] := z'(£j) 
with fj = jk and k being the step size, we have the following quadrature approximation 


J CXJ 

( 4 . 5 ) U h (t) = — e^cpjZjVh, 

J = - oo 

and the truncated quadrature approximation 

k N 

( 4 - 6 ) U Nih (t) = — eZi % z 'j, 

ni j=-N 

with cpj = (zjiv(zj)I + Ah) 1 w(zj)vh■ To compute UN,h(t), we need to solve only 7V + 1 elliptic problems, 
instead of 2N + 1 elliptic problems, by exploiting the conjugacy relations: Z-j = ~z], w(z-j) = w(z-j), 

(p_j = cpj, j = 1, • • • , N. Indeed, since z'j = z'(£j) = iAcos(i^ — ip), denoting by Q = Acos(i — ip), (4.6) 
is reduced to 


(4.7) 


UN,h{t) = — 


N 


-e^oCo 

3 =1 


Hence we solve the following complex-valued elliptic problems 


(4.8) (zjw(zj)I + A h ) <pj = w(zj)v h , j = 0 , ...,1V. 

These problems are independent of each other and can be solved in parallel, if desired. 
Next, we define a strip <S a> t, C C by 


S a ,b — {p = £ + ii? : for all £ € R and y £ (— 6 , a)}. 


The following lemma recalls a known error estimate for the quadrature [29] [44, Theorem 2.1]. The 
quadrature is exponentially convergent, provided that the integrand g is analytic on a strip S a ,b with 
some additional conditions. 
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Lemma 4.1. Let g be an analytic function in a strip S a ^ for some a, b > 0, and I and Ik, for k > 0, be 
defined by 

/ OO 00 

g{x) dx and I k = k ^ g{jk), 

respectively. Furthermore, assume that g{z) —> 0 uniformly as \z\ —> oo in the strip S a and that there 
exist M+ > 0 and M_ > 0, which may depend on a and b such that 

/ OO /-OO 

\g{x + ir)| dx < M + , lim / \g{x — \s)\dx < M_. 

-oo s ~>- b ~ J — oo 

Then the approximation error can be bounded by 

\I-Ik\ <E++E~, 

where 

E + = — 7 — 77 ~—- and E~ = 




g ‘la'KI k ^ g2bn/k ^ 

The next lemma gives one crucial estimate on the map z(jp) over the strip S a .b- Even though the 
hyperbolic contour (4.1) has been extensively used, the estimate on the map z{p) below seems to be new 
and it is of independent interest. 

Lemma 4.2. Let p = £ + ip with £, g £ 1R. Then with a = 7 r /2 — ip — e and b = if — e, for small e > 0, 
there holds 

z'{p) 


(4.9) 

(4.10) 


z{p) £ and 

z{p) £ and 


z{p) 

z'{p) 


*(P) 


< 


< c 


Vp £ 5a,0; 

Vp £ So,b- 


Proof. For p = £ + ip with f, g £ K, then the image z{p) in the parameterization (4.1) is given by 
z{p) = A(1 — sin(^ + g) cosh(£)) + iA cos {ip + g) sinh(£), 
and its derivative z'{p) is given by 

z'{p) = Acosh^cos(V’ + g) — i sinh £ sin(i/> -f g). 

By writing z = x + \y, it can be expressed as the left branch of the hyperbola 

x - A \ l y 


A sin^ + g) 


cos{ip + g) 


= 1 . 


It intersects the real axis at x = A(1 — sin{ip + g)) and has the asymptotes y = ±{x — A) cot(^ + g). Next 
we show the estimates (4.9) and (4.10). First, for p £ 5 0j o, i.e., g £ [0,a], z{p) lies in the sector 
Using the elementary identity sinn x = cosh 2 x — 1, the fact <p := g + if £ {ip, tt/2 — e), and the estimate 
sin(7r/2 — e) ~ 1 — e 2 /2 < 1 — e 2 /3 for small e, we have for all £ £ R 

2 


cosh 2 (£) — sin 2 (<p) 


z'{p) 

2 

cos(<p) cosh(£) — isin(p) sinh(£) 

z{p) 


(1 — cosh(£) sin(p)) + isinh(£) cos(tp) 


cos 2 (<p) cosh 2 (£) + sin 2 (p) sinh 2 (£) 


1 — 2 cosh(£) sin(<p) + cosh 2 (£) sin 2 (p) + sinh 2 (£) cos 2 (p) (cosh(^) — sin(p )) 2 


cosh(£) + sin(p) 1 + sin(<p) 


6 

< —77 ■ 


cosh(£) — sin(p) 1 — sin(<p) 1 — (1 — e 2 /3) e 2 

Hence the estimate (4.9) holds true. Now we turn to the case p £ S o,&, i.e., g £ [—b, 0]. Then z{p) lies in 
the sector S 7r _( I)+ ^) C Y K _ e . Further, by noting := g + ip £ (e, ip), we have for all £ £ R. 


z'{p) 


z{p) 


< 1 + sin(<p) 1 + sin(^) 
— 1 — sin(p) — 1 — sin {ip) 


Then the desired result (4.10) follows directly. 


□ 
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The next result gives the analyticity of and an estimate on the integrand g(£,t) on the strip S a ,b- 

Lemma 4.3. Let p = £ + i ?7 with £, 77 £ M and g{p,t) be defined by (4.4). Then fjip,t) is analytic on the 
strip S a ,b, and the following estimate holds: 

\\g(pM < ^e A ( 1 - sin ^) cosh «)) t ||^|U 2 ( o) Vp G s a , b . 

Proof. For p = £ + ip with £, 77 G M, the image z(p) in (4.1) is given by 

z{p) = A(1 — sin(^ + rf) cosh(£)) + iA cos{ip + 77 ) sinh(£). 

By Lemmas 4.2 and 2.1, and Remark 2.1, z(p)w{z(p)) G S OT _ e /, with e' > 0. Hence the function 

g(p,t) = e z ( p)t {z(p)w{z(p))I + A^ 1 w(z(p))z'(p)v h 

is analytic in the strip S a ,b- It remains to show the estimate. First, we consider the case p £ So,b- By 

(4.10) , z{p) £ Stt-e- Then, by Lemma 2.1 and Remark 2.1, z(p)w(z(p)) G E,r_ e /, with e' = ce. By the 
resolvent estimate ( 2 . 2 ), we deduce that for small e > 0 , there holds 

(4.11) || {zl + Ah) 1 || < c/|3(z)| < c/| 2 sin( 7 r — e)| < c/{\z\e') Vz G S( r _ e ,. 

Meanwhile, for any p G S o,&, there holds 

3?(z(p)) = A(1 - sm(ip + 77 ) cosh(£)), 

which together with the resolvent estimate (4.11) and Lemma 2.1 yields 

\\g(p,t)\\ < ce st( - z{p))t \z , (p)w{z(p))\ \\(z(p)w(z(p)) + A)" 1 ! 

Ap) I 


IMU 2 (fi) 


< -e 
e 


£ „A(l-sincosh(£))£ 


zip) 


IKIU 2 (Q)- 


This together with (4.10) yields the desired assertion. The case p £ <S 0j o is more direct. Then (4.10) 
and Lemma 2.1 imply that zip)w{z{p)) G Eg/ with 9' G (7 t/2 , 7r) depending only on ip. Then the desired 
assertion follows from (4.9) and the resolvent estimate (2.2). □ 

Now we can give an error estimate for the quadrature approximation 

Lemma 4.4. Let Uh{t) and UN,h{t) be defined in (4.3) and (4.6), respectively, and the contour be para¬ 
metrically represented by (4.1). Then with the choice k = Co/N and A = ciN/t, there holds 

II Uh(t) - t/jv,/i(t)||.L 2 (n) < ce~ c ^IMIl^o), 
where the constant c and c' depend on the choice of ip in (4.1). 

Proof. We use the following splitting 

Uh ~ Uw,h = {uh — Uh) + (Uh — UN,h) ='■ E q + E t , 

where E q and E t denote the quadrature and truncation error, respectively. We apply Lemma 4.1 to bound 
\\E q 11/, 2 (q) - To this end, we set a = 7 r /2 — ip — e and b = ip — e. For p = £ + ia, zw(z) lies in the sector 
Eg for some 9 G ( 7 t/ 2 , 7 r). Note the elementary inequalities cosh £ > 1 + £ 2 /2 and 1 — sin( 7 r /2 — e) < e for 
small e > 0. These together with the choice A = c\N/t and Lemma 4.3 yield 


r 


|g(£ + ia)|d£ 


l 2 ( n) 


c f°° 
<~ e 
e Jo 


ci N(l— sin(7r/2— e) cosh(£)) 


^||^hlU 2 (a) 


< -e 
e 


c\N e 


f 


0 —c\N sin(7r/2—e)£ 2 /2 


d£ph\\ 


L 2 {Q) 


< -N-h^WvhWwy 


Using Lemma 4.1, for k = cq/N we have 


\E+\\ L 2 [n) < -IV-^e-l 277 ^/ 2 -^ 6 )/ 00 " 60 !)^. 
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Next we bound the error due to the lower half. For the choice p = £ — \b, A = ciN/t and appealing again 
to the inequality cosh£ > 1 + £ 2 /2, we deduce 


|<7(£-i6)|d£ 


L 2 (fi) 


< - / e ^(i-“Wcodi(0)^ ||t, h || ia(n) 

6 Jo 

/»oo 

< / e -ciiVsin( e )« 2 /2^ || Wh || ia(n) 

6 Jo 


c 

^3/2 J 

Then for the choice k = cq/N, Lemma 4.1 yields the following estimate 

||£-|U 2 (q) < 

Further, by using cosh(£) > cosh(co) + sinh(co)(£ — cq) for £ > cq, the truncation error ||£7 t ||z, 2 ( 0 ) can be 
simply estimated by 


r r°° 

\E t \\ L 2 {Q) < - / gCiAr(l-sin(^) cosh({)) ^ || 

6 j Co 

r I°° 

< _ e ciA(l—si n (V')cosh(co)) / g — ciIVsin(i/>) sinh(co)(f —Co) 

“ 6 J cn 


d£ll v /tlU 2 (n) 


-iV- 1 e ClJV[1 " sinWcosh(ct,)1 ||i; /l || L 2 (n) . 


Ci sin(-0) sinh(co)e 

Finally, by disregarding e terms, balancing asymptotically the exponential parts in || E+ ||x, 2 (f 2 ), \\Eq ||L 2 (n) 
and \\E^ ||x, 2 (fj), we arrive at 

2tt(tt/2 — i/j)/co = 2-Kip/cQ — Ci = —ci(l — sin(^) cosh(co)). 

We may express the parameters cq and ci in terms of ip: 


cq = cosh 


27 rip 


and ci = (47 Tip — tt 2 )/ cosh 


27 Tip 


(Anip — 7 r 2 ) sin ip J 


k (47 Tip — 7T 2 ) sin^i. 

Finally we minimize the ratio 

B{ip) = Ci- 2mp/c 0 

with respect to the parameter ip, which achieves the minimum at ip = 1.1721 and hence, 

c 0 = 1.0818, Ci = 4.4920 and B{ip) = -2.32, 

which are identical to those values given in [44] . Then collecting the balanced asymptotic bound and the 
rest from \\E+\\ L 2 {n) , ||-E“|| L 2 (n) and ||^ t "|U 2 (fi) yields 

K(t) - U N , h (t)\\ L2{n) < c (e - 1 ^- 1 / 2 + e-^N- 1 ' 2 + e-^ 2 ~^^ +c ^ N \\v h \\ L 2 {n) . 

Now by choosing e = 1/iV, we get 

II Uh{t) - UN,h(t)\\L 2 (u) < ce ( '~ 2 32+ ™ )w ||w/i||i 2 (n)- 

which together with the fact (logx)/a; < 1/e for x > 1 and the Testability of the projection Ph yields 
the desired result. □ 

Last, we give the main result of this section, i.e., error estimates for the fully discrete scheme (4.6). It 
follows from Theorems 3.1 and 3.2, and Lemma 4.4 and the triangle inequality. 

Theorem 4.1. Let u(t) be the solution of problem (1.1), and Ujsi,h(t ) be the quadrature approximation 
defined in (4.6), with the parameters chosen as in Lemma f-4- Then with l\ (t) = (log2T/f) -1 , the 
following estimates hold. 

(a) If Av € L 2 (Li) and Vh = RhV, then 


II u{t) - U Nth (t)\\ L 2 {n) < c (e C ' N + h 2 ^ \\Av\\ L 2 (n) . 
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(b) If v € L 2 {H) and Vh = PhV, then 

IK*) - U Nth ,\\ L *(si) < ct (e~ c ' N + ||w|| L 2(n). 


5. FULLY DISCRETE SCHEME II: CONVOLUTION QUADRATURE 

Now we develop a second fully discrete scheme based on convolution quadrature generated by the 
backward Euler method, and show that the scheme is first order convergent. 


5.1. Time stepping based on convolution quadratures. To describe the fully discrete scheme, we 
divide the interval [0,T] into a uniform grid with a time step size r = T/N , N € N, with 0 = to < *i < 

... < tjv = T, and t n = nr, n = 0,..., N. The general construction of convolution quadrature is as 
follows [25, 7]. Let (a,p) be a stable and consistent implicit linear multistep method, with (ct, p) being 
its characteristic polynomials. Then we define convolution quadrature weights {bj}fjL 0 by the expansion 
coefficients of 

=({)=£«'=/ (^f) (•(«>*»• 

We consider only the simplest case, i.e., the backward Euler method, for which the convolution quadrature 
weights {bj}f i 0 are defined by 

(5.1) S«) = £>,?=/‘(L^ 

U V T 

The convolution quadrature weights {bj} can be computed efficiently using the fast Fourier transform [37], 
in view of Cauchy’s theorem. Then the convolution quadrature Q T p for a Riemann-Liouville fractional 
derivative R Dfp := f,j(t — s)~ a ip(s)ds generated by the backward Euler method is given by 

n 

(5.2) (Qr<p)(tn) = b n -jip(jT). 

3 =0 

Following this general construction, we now derive the time stepping scheme. The approximation 
Qn(<p) to the Riemann-Liouville fractional derivative R D{*(p(t n ) is given by [7, 16]: for any n = 1, 2 ,..., N: 

n 

(5.3) Qn(p) = ^ bn—jPftj ) ; 

3 =1 

where the weights {bj} are generated by (5.1). Recall also the defining relation of the Caputo derivative 
using the Riemann-Liouville derivative [20, p. 91, equation (2.4.4)] Dfu = R Df{u — u( 0)). Upon applying 
the convolution quadrature to the term on the right hand side and using it for the semidiscrete problem 
(3.2), we arrive at the following fully discrete scheme for the model (1.1): for n = 1, 2,..., TV 

(5.4) Q n (U h ) + A h U% = Q n (l)v h , 

with £7° = Vh- Throughout, we denote the generating function /3 of a sequence {/3j}^L 0 by /?(£) = 

ofot 3 - 

Remark 5.1. Compared with the general construction (5.2), the term corresponding to j = 0 is omitted 
in our fully discrete scheme (5.4). This choice was taken earlier in [26, 3]. 


p(a) da = 


i-£ 


w 


i-£ 


5.2. Error analysis. Now we carry out the error analysis of the fully discrete scheme (5.4), following 
the strategy outlined in the pioneering work [26]. To derive L 2 (fl)-error estimates, we split the error into 

e” = u(t n ) - Uf; = ( u(t n ) - u h (t n )) + (uh(t n ) - Uh). 

In view of Theorems 3.1 and 3.2, it suffices to establish a bound on | \uh(t n ) — UJ}\\ L 2 ^ n y The proof relies 
on the following splitting 

Uhitn) - K = Vh(t) - Y£, 







ERROR ESTIMATES FOR DISTRIBUTED ORDER FPDE 


15 


where 

Vh(t ) = u h {t) - v h and Yjf = - v h . 

First, we derive representations of the semidiscrete solution yh and fully discrete solution Y ) t . 
Lemma 5.1. Let the kernel K(z) be defined by 

(5.5) I<(z) = -z _1 ( 2 w( 2 :)/ + A^^Ah 

and x( z ) = —• Then yh and Yff can be represented by 

Vh(t) = 2^1 / e zt K(z)v h dz and Y h = YJ f e ztr '- 1 K{x{z))v h dz, 

respectively, with the contour V T = {z G T#^ : |3(z)l < 7r/r}. 

Proof. By its definition, yh satisfies the problem: 

d[ m1 j )h + A h y h = -A h v h , 
with yh(0) = 0. The Laplace transform gives 

zw(z)y h (z) + A h y h (z) = -z~ x A h Vh- 

Hence, yh{z) = K(z)vh, with K(z) = —z~ 1 (zw(z)I + Ah) -1 Ah, and the desired representation for yh(t) 
follows from the inverse Laplace transform. Next, the fully discrete solution Yff satisfies the following 
time stepping scheme 

Qn(Xh) + AYff = —A h Vh, 

with Y° = 0. Now multiplying both sides by summing from 1 to oo and noting Y® = 0 yield 

OO 

Qn(Y h )C + A h Y h (0 = -£/(l - OA h v h . 

n— 1 

Using the condition Y® = 0, we have 

oo oo n 

E Qn(Y h )C = E E ( b n-iC~ j ) (y£?) = ((1 - 0/r)w(( 1 - 0/r)Y h (0- 

n— 1 n— 0 j= 0 

Thus, by simple calculation, we deduce 

Y h (0 = (£/T)K(( 1-O/tK, 

and it is analytic at £ = 0. Then Cauchy theorem implies that for g small enough, there holds 

Y h n = ^-[ C n K{{l-0/T)v h dti. 

2rm J , € | =e 

Now, by changing variable £ = e~ ZT , we obtain 

Y h = J T e zt ^K((l^e~n/r)v h dz, 

where the contour T 0 = {z = — In {q)/t + i y : |y| < 7r/r} is oriented counterclockwise. We obtain the 
desired representation by deforming the contour To to T r = {z G Tg^ : |S(z)| < 7r/r} and using the 
periodicity of the exponential function. □ 

By Lemma 5.1, we can write the difference between Yff and yh(t n ) as 

Vh(t n )-Y£ = I + II, 

where the terms / and II are given by 

(5.6) / = [ e ztn K(z)v h dz 
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and 

(5.7) II = f e 2t " (K(z) - e~^K(x(z))) v h dz. 

This splitting is essential for the error analysis below. Since the function |e _2T | is uniformly bounded on 
the contour T r , we have 

II K(z) - e~ ZT K(x(z))\\ < \e~' T \\\K{z) - K( X (z ))|| + |1 - e— \\\K{z)\\ 

(5.8) < c\\K(z) - K{ X (z))\\ +ct\z\\\K(z)\\ 

< c\\K(z ) - K(x{z))\\ + ct , 

where the last line, using the resolvent estimate (2.2), follows from the inequality 

n^)ii = kri - /+ zw(z)(zw(z )+^)- i n < ci^- 1 . 

Hence, it remains to bound the term ||A'( 2 ) — K(xi z )) || t which will be carried out in several steps. 
First we recall a bound on the function X (z) = r _1 (l — e~ ZT ) [17, Lemma 3.1]. 

Lemma 5.2. Let X (z) = r _1 ( 1 — e~ ZT ). Then for all z £ T r , there hold 

Ix(-) ^ ~l < c\z\ 2 t and Cx\z\ < \ X (z)\ < c 2 \z\, 
and x(z) lies in a sector Eg/ for some 6' £ (tt/2,tt). 

Next we give one crucial error estimate on the approximation x(z)w(x(z)) to the kernel zw(z). 
Lemma 5.3. For z £T t , the following bound holds: 

\x{z)w{x{z)) - zw{z )| < ct\z\ 2 w(\z\). 


Proof. By the intermediate value theorem, for z £ T t , we have 

fX(z) 

\ X (z) a - z a \ = a 


„a-l 


ds 


< a\x(z) — z | max | z., 


|a-l 


r)e[o,i] 


where z v = rjx{z) + (1 — rf)z with rj £ [0,1]. Next we claim | z. 


77 1 1 < c\z 


-l 


for 77 £ [0,1]. To this end, we 


split T r into T r = T+ U TJ; U T t , with T^ being the rays in the upper and lower half plane, respectively, 
and T% is the circular arc. For z £ F°, by the Taylor expansion of e~ ZT , we have 

/ OO 

z v = z\ 1 + ? 7 ^(-l) J 


j =1 


U + 1 )! 


In view of the trivial inequality |zt| < 1 for z £ we deduce l^l -1 < c| 2| _1 for z £ F°. It remains to 
show the assertion for z £ T+, and the case z £ T~ follows analogously. First we show Q(x(z)) > 0 for 
z £ r+. For z = re‘^ K ~ e l with rr £ (<5, 7 r/ sin0) we have 


x( z ) = - (1 - 


g rr cos 0 g—irr sin 


in 0^ 


and therefore using rrsin0 < 7 r, we get Sy(x(^)) > 0. Then Lemma 5.2 yields 


Q 

\z v | > min(|z|, |x(z)|)cos - > c|z|. 


This shows the desired claim. Hence, appealing to Lemma 5.2 again implies that for z £T t there holds 
(X(z) a - z a )fx(a) da 

which concludes the proof of the lemma. □ 


f 


< f \x{z) a — z a \n(oi) da < ct\z\ f \z\ a ^(a) da = ct\z\ 2 w(\z\), 

Jo Jo 


Next we give a crucial error estimate on the approximation K(xi z )) to the kernel function A'(z). 
Lemma 5.4. Let x( z ) = (1 — e~ ZT )/r. Then for the kernel K{z) in (5.5), there holds 

\\K(z) - K(x(z))\\ < ct Vz £ r r . 
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Proof. Let B(z) = zK[z). Simple computation shows 

B (z) - B (x(z)) = zw(z) (zw(z)I + A^ 1 - x(z)w(x(z)) (x(z)w(x(z))I + A h ) _1 
= zw(z) (( zw(z)I + A h y x - (x{z)w{x(z))I + A h )~^ 

+ (zw(z) - x(z)w(x(z))) {x(z)w{x(z))I + Ah) -1 := I + II. 

First, by Lemmas 2.3 and 5.2, there holds 

\x{z)w{x{z))\ > c|x(V)Klx(>)|) > c\z\w{\z\). 

Further, by Lemma 2.1 and (2.2) and Lemma 2.3, we have 

(5.9) || (zw(z)I + Ah)- 1 1| < c\zw(z)\~ l < c(| 2 |u;(| 2 ;|)) _1 . 

Likewise, in view of Lemmas 5.2 and 2.1 and (2.2), we have 

(5.10) \\(x(z)w(x(z))I + Ah)- 1 1| < c\ x (z)w(x(z))\- 1 < c(|^|w(|^|)) _1 . 

Now, by the identity 

(zw{z)I + Ah)' 1 - (x(z)w(x{z))I + Ah) -1 
= ( zw(z ) - xO*MxO*))) (zw{z)I + A h y l (x{z)w(x{z))I + Any 1 , 

Lemma 2.3, (5.9) and (5.10), the first term I can be bounded by 

||/|| < ct\z\ 3 w(\z\) 2 \\(zw(z)I + A h y l \\\\(x{z)w{x(z))I + Ay-^l < ct\z\. 

Likewise, with Lemma 5.3 and (5.10), the second term II can be bounded by 
||//|| < \zw{z) -x{z)w{x{z))\\\{x{z)w{x{z))I + A h y l \\ 

< cr 1 | 2 xt;( 1 |)|(|-s |)| —1 < cr| 2 :|. 

Hence we bound || B(z) — B{x{z)) || by 

\\B(z)-B( X (z))\\<ct\z\. 

Last, by Lemma 5.2 and ||H(z)|| < c, we bound || K(z) — K(xi z ))\\ by 

II K(z) - K( X (z)) || < l^" 1 - xizy^WB^W + W-'WBiz) - B( X (z))\\ 

< c\z - x{z)\\z\~' 2 + CT < CT, 

which completes the proof of the lemma. □ 

Now we can state an error estimate on the time discretization error for nonsmooth initial data, i.e., 
v e L 2 yi). 

Theorem 5.1. Letuh andUfi be the solutions of problems (3.2) and (5.4) withv € L 2 yi), U° = Vh = PhV 
and f = 0, respectively. Then there holds 

\\u h {t n ) - K\\ L 2 {Q) < cr^ 1 ||u|| L 2 (n) . 


Proof. It suffices to bound the terms I and II defined in (5.6) and (5.7), respectively. We choose 5 = ty 
in the contour Tye. By (2.2) and direct calculation, we bound the first term I by 


(5.11) 


pOO 

L 2 {Q) <c e rt n co S 9 r -l dr\\v h \\ L 2(n) 

J TV/(t sin 6) 


< CT 


L 2 (Q) 


f 


e rt n cos 9 dr < CTt -l 


L 2 (Q.)- 


Using Lemma 5.4, we arrive at the following bound for the second term II : 


(5.12) ||7/||l2(q) < ct||u/j||l2(q) (j 


7r/(r sin 0) 


' dr + J e^tydA < ct n 1 r||u? l || L 2(Q). 
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Combining estimates (5.11) and (5.12) yields 

II Vh{tn) - Yh llz, 2 (n) < CTtn 1 II II L 2 (n), 

and the desired result follows directly from the identity U% — Uh(t n ) = Yjf — yh(t n ) and the stability of 
the projection Ph in L 2 (H). □ 

Remark 5.2. The L 2 (Ll) stability of the time stepping scheme (5.4) follows directly from Theorem 5.1. 

Next we turn to smooth initial data, i.e ., Av £ L 2 (Ll). To this end, we first state an alternative 
estimate on the solution kernel K(z). 

Lemma 5.5. Let K s (z) = — z~ 1 (zw(z)I + Ah)'. Then for any z £ T t , there holds 

\\K'{z)-K'(x(z))\\<ct^-. 

\ z \ - 1 

Proof. Let B s (z ) = —(zw(z)I + A^) -1 . Then by the trivial inequality 

B s (z) - B s ( x (z)) = X (z)w( X (z)) - zw(z)) ( zw(z)I + A ,)” 1 ( X (z)w( X (z))I + Ah)- 1 
Lemma 5.3, and (5.9) and (5.10), we deduce immediately 

\\B s {z)-B s ( X m\<cT\w{z)\- 1 . 

Appealing to 5.9 again, we have ||.B s (.g)|| < c)zw(z)\~ 1 , and thus 

II K s (z) K°( X {z)) II < |~ _1 ^ X(*r 1 |||£ s (z)ll + IxizT^iz) - B°(x(z))\\ 

< c\z — x(^)||~| — 3 |^u (^)| —1 + ct\zw{z)\^ 1 < ct| 2 u;(.z)| _1 . 

Then the desired result follows from Lemma 2.3. □ 

Now we can state an error estimate for smooth initial data Av £ L 2 (Ll). 

Theorem 5.2. Let Uh and UJ) be the solutions of problems (3.2) and (5.4) with Av £ L 2 (Q), U® = Vh = 
RhV and f = 0, respectively. Then for ^(t) = log (max(f - 1 ,2)), there holds 

II u h {t n ) - U% || L 2 (n) < cT£ 2 (t)\\Av\\ L 2 {n) . 


Proof. Let K s (z ) = 1 (zw(z)I + Ah) . Then we can rewrite the error as 

Vh(tn) ~Yff =-^~-f e ztn K s (z)A h v h dz 

2lT1 J r«, 4 \rr 

+ J e ztn ( K\z) - e~ ZT K s ( X (z ))) A h v h dz := I + II. 


(5.13) 


By Lemma 5.5 we have for ;6 T t 


II K s (z)-e ZT K s (x{z))\\ < ct 


log I ^ I 
\z\-l 


By setting <5 = 1 jt n and by the monotonicity of the function f(x) = ^ on R + , we derive the following 

bound for the term II 

ptc /(r sin 6) 


II^IU 2 (n) < cT\\A h Vh\\L 2 {Q) 


jrt n cos 6 


1 /tn 


logr 
r — 1 


dr+ f e c°8^ 1 °g(*n 1 ) d A < c 
J —6 1 tn J 


lQg(*n ) 

1-tn 


T\\A h Vh\\ 


L 2 


Now (2.2) implies that for all 2 £ Tg^, ||i^ s ( 2 ;)|| < c|^| 1 \zw(z)\ 1 . Therefore, using Lemma 2.3, we 
deduce 


(5.14) 


POO 

U 2 (n) < c\\A h v h \\ L *{n) / e rt " cose r' 2 |u;(r )|^ 1 dr 

J 7 r/(r £ 


t/(t sin 6) 


POO 

< CT\\A h v h \\ L 2 {n) / e 

Jl/t„. 


l0g r ■ dr < A 0 ^ 1 ). 


r — 1 


1 -t„ 


\\AhVh\\L 2 (n)- 
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Finally, we observe that if t n x > 2, i.e. t n < 1/2, then lo ^ t TI 1 < 21og(t n 1 ). Otherwise if t rl 1 < 2, i.e. 

t n > 1/2, then by the monotonicity of the function f(x) = on R + , we deduce lo ^” ^ < 

21og(2). Then the desired result follows from (5.2), (5.14) and the identities Ufi — Uh{t n ) = Yff — yh{t n ) 
and A h R h = P h A. □ 

The next theorem gives error estimates for the fully discrete scheme (5.4), which follow from Theorems 
3.1, 3.2, 5.1 and 5.2 and the triangle inequality. 

Theorem 5.3. Let u and be the solutions of problems (1.1) and (5.4) with U® = Vh and f = 0, 
respectively. Then for t\(t) = log(2T/t) _1 and I 2 (t) = log (max(t -1 ,2)) and t n = nr, the following error 
estimates hold. 

(a) If Av £ L 2 ffV) and Vh = RhV, then for n > 1 

\\u(t n ) - UJ/\\ L2(n) < c(rt 2 (t n ) + h 2 )\\Av || L 2 (n) . 

(b) If v £ L 2 ffl) and Vh = PhV, then for n > 1 

I \u(t n ) - U£ || L 2 (n) < c T (t + h 2 h(t n )) t- 1 ||u|| I/ 2 (n) . 

Remark 5.3. For distributed order time fractional diffusion, the error estimate involves a log factor in 
time for smooth initial data, which is reminiscent of the asymptotic behavior of the solution at small time, 
cf. Theorem 2.1 . This factor is not present for the single term and multi-term time fractional diffusion 
[14, 15], 


6. Numerical experiments and discussions 

Now we present numerical results to verify the convergence theory. To this end, we let the domain Q 
to be the unit internal fl = ( 0 , 1 ) and consider the following three examples with smooth, discontinuous, 
and singular initial data: 

(a) v(x) = sin(27rx) £ H 2 (fl) n Hq (fl); 

(b) v = X(o,i/ 2 ) G H 1 / 2 ~ e (fl) with e £ (0,1/2), and xs the characteristic function of a set S; 

(c) v(x) = a : -1 / 4 £ H 1 / 4 ~ e (Tt) with e £ (0,1/4). 

We measure the temporal discretization error by the normalized L 2 ffX) errors \\u(t n )— tlN,/i(tri)||z, 2 (n)/IMlL 2 (n) 
or || u(t n ) — Cl^ I ||i, 2 (n)/||i ; ||L 2 (n) J and the spatial discretization error by the normalized L 2 (ffL) and H 1 ( fl) 
errors, i.e., )|u(t) - u h (t)\\ L 2 ^/\\v\\ L 2 ^ and ||V(w(f) - u h (t))\\ L 2 ^/\\v\\ L 2 ^ n y In the computations, we 
divide the domain H into M equally spaced subintervals with a mesh size h = 1/M. Since the exact 
solution u{t) is not available in closed form, we compute the reference solution using a much finer mesh. 

6.1. Numerical results for the semidiscrete scheme. First we examine the convergence behavior 
of the space semidiscrete scheme. To this end, we fix N = 10 in the Laplace transform approach such 
that the error due to time discretization is negligible. The numerical results are given in Tables 1-3. In 
the table, rate denotes the empirical convergence rates when the mesh size h halves, and the numbers 
in the bracket denote the theoretical rates. For all three initial data, the L 2 {Ll) and H 1 (Cl) norms of the 
error exhibit second and first order convergence rates, respectively, which agrees well with the theoretical 
prediction, cf. Theorems 3.1 and 3.2. The convergence of the semidiscrete scheme is robust in that the 
convergence rates hold for both smooth and nonsmooth initial data. The error increases as t —> 0, which 
is attributed to the weak singularity of the solution as t —> 0, cf. Theorem 2.1. 

6.2. Numerical results for the fully discrete scheme I. Next we illustrate the convergence of the 
first fully discrete scheme based on the Laplace transform. To make the spatial discretization error 
negligible, we fix the spatial mesh size h at h = 10~ 5 . In all numerical simulations, the optimal contour 
parameters A and if in the parameterization (4.1) and k in (4.7) are chosen as suggested in the proof 
of Lemma 4.4 (see also [44]). Moreover, A is fixed, independent of t, with which the elliptic problems 
(4.8) are solved for each time t. The numerical results are summarized in Tables 4 and 5 for the weight 
functions p(a) = (a — 1/2) 2 and /r(a) = X[i/ 2 ,i](<Y), respectively. The results indicate an exponential 
convergence with respect to the number N of quadrature points on hyperbolic contour, decaying at 
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Table 1. Numerical results for the standard semidiscrete Galerkin FEM for smooth 
initial data, Example (a) with N = 10 and ^t(a) = (a — 1/2) 2 . 


t 

M 

10 

20 

40 

80 

160 

320 

rate 

1 

L 2 (n) 

2.79e-5 

7.02e-6 

1.76e-6 

4.39e-7 

1.09e-7 

2.70e-8 

2.00 

( 2 . 00 ) 


8.84e-4 

4.44e-4 

2.22e-4 

l.lle-4 

5.23e-5 

2.36e-5 

1.05 

( 1 . 00 ) 

10- 1 

L~(VL) 

2.40e-4 

6.05e-5 

1.52e-5 

3.79e-6 

9.35e-6 

2.33e-7 

2.00 

( 2 . 00 ) 

H L \n) 

7.03e-3 

3.53e-3 

1.77e-3 

8.84e-4 

4.16e-4 

1.88e-4 

1.00 

( 1 . 00 ) 

10 -3 

L 2 (n) 

6.38e-3 

1.61e-3 

4.03e-4 

1.01e-4 

2.51e-5 

6 . 21 e -6 

2.00 

( 2 . 00 ) 

i njTf 

1.41e-l 

7.04e-2 

3.53e-2 

1.76e-2 

3.75e-3 

1.65e-3 

1.07 

( 1 . 00 ) 


Table 2. Numerical results for the standard semidiscrete Galerkin FEM for discontin¬ 
uous initial data, Example (b) with N = 10 and ^t(a) = (a — 1/2) 2 . 


t 

M 

10 

20 

40 

80 

160 

320 

rate 

1 

L 2 (n) 

3.97e-5 

9.94e-6 

2.48e-6 

6.21e-7 

1.55e-7 

3.87e-8 

2.00 

( 2 . 00 ) 

1 T T JTiy 

1.26e-3 

6.29e-4 

3.15e-4 

1.55e-4 

7.63e-5 

3.68e-5 

1.01 

( 1 . 00 ) 

10" 2 

T 2 ( 0 ) 

5.81e-4 

1.45e-4 

3.64e-5 

9.12e-6 

2.28e-6 

5.69e-7 

2.00 

( 2 . 00 ) 


1.28e-2 

6.38e-3 

3.19e-3 

1.57e-3 

7.73e-4 

3.73e-4 

1.02 

( 1 . 00 ) 

10 -3 

L~(VL) 

6.34e-3 

1.59e-3 

3.96e-4 

9.92e-5 

2.48e-5 

6.18e-6 

2.00 

( 2 . 00 ) 

H l {n) 

1.73e-l 

8.65e-2 

4.32e-2 

2.14e-2 

1.04e-2 

5.06e-3 

1.02 

( 1 . 00 ) 


Table 3. Numerical results for the standard semidiscrete Galerkin FEM for singular 
initial data, Example (c) with N = 10 and //( a ) = (a — 1/2) 2 . 


t 

M 

10 

20 

40 

80 

160 

320 

rate 

1 

L' 2 (Q) 

3.82e-5 

9.67e-6 

2.44e-6 

6.12e-7 

1.53e-7 

3.79e-8 

2.00 

( 2 . 00 ) 

H l {n) 

1.21e-3 

6.13e-4 

3.09e-4 

1.55e-4 

7.33e-5 

3.33e-5 

1.05 

( 1 . 00 ) 

10" 2 

L 2 (n) 

6.72e-4 

1.69e-4 

4.23e-5 

1.06e-5 

2.63e-6 

6.51e-7 

2.00 

( 2 . 00 ) 


1.38e-2 

6.92e-3 

3.47e-3 

1.74e-3 

8.18e-4 

3.71e-4 

1.06 

( 1 . 00 ) 

10“ 3 

L 2 (n) 

3.48e-3 

8.76e-4 

2.20e-4 

5.49e-5 

1.37e-5 

3.36e-6 

2.00 

( 2 . 00 ) 

H^n) 

1.49e-l 

7.45e-2 

3.73e-2 

1 . 86 e -2 

8.76e-3 

3.97e-3 

1.07 

( 1 . 00 ) 


a rate about e ~ 2 ' 15N and e _214Ar for n{a) = {a — 1/2 ) 2 and n(a) = X[i/ 2 ,il ( a ), respectively, which 
agree well with the theoretical predictions from Theorem 4.1. Note that even though the weight function 
p(a) = Xri/ 2 , 1 ] i a ) does not satisfy the assumption p(0)^(l) > 0, the empirical convergence rates still agree 
well with the theoretical prediction, which calls for further theoretical study. Further, the convergence 
rate is independent of time t, and thus the scheme is robust. Interestingly, the smoothness of the initial 
data v does not affect much the time discretization errors, even for small time instances, cf. Table 6 . 

One salient feature of the fully discrete scheme I is that it allows computing the solution at any 
arbitrarily large time directly. This allows one to examine the asymptotic behavior of the solution as the 
time t — > oo; see Table 7 and Fig. 1. In particular, one clearly observes the logarithmic decay of the 
solution, as predicted by [22, Theorem 2.1]; see also Fig. 1. This numerically verifies the ultraslow decay 
asymptotics for distributed order diffusion process, in comparison with sublinear decay for subdiffusion 
and exponential decay for normal diffusion. 

6.3. Numerical results for the fully discrete scheme II. Last we verify the convergence of the fully 
discrete scheme II, i.e., convolution quadrature based on the backward Euler method. By Theorem 5.3, 
it exhibits a first order convergence with respect to the time step size r. This is fully confirmed by the 
numerical results in Tables 8 and 9 for the weight functions n{a) = (a — 1/2 ) 2 and /J.(a) = X[i/ 2 ,il(c^)? 
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Table 4. The L 2 errors for initial data (a)-(c) with h = 10” 5 and n{a) = (a — 1/2) 2 , 
by the Laplace transform method. The notation r denotes the exponential convergence 
rate in the error || u% >h - u(t n )\\ L 2 ^ < Ce~ rN . 


case 

t \ N 

3 

5 

7 

9 

11 

13 

r 


1 

1.33e-6 

1.49e-8 

1.26e-10 

2 .20e-12 

3.54e-14 

8.24e-17 

2.35 

(a) 

10" 2 

4.78e-6 

7.36e-7 

2.77e-9 

5.45e-ll 

4.88e-13 

2.23e-14 

1.92 


10“ 3 

8.30e-5 

8.78e-7 

3.81e-9 

7.55e-ll 

6.43e-13 

1.23e-14 

2.26 


1 

3.34e-6 

3.56e-8 

2.85e-10 

5.76e-12 

8.68e-14 

1.25e-15 

2.17 

(b) 

10~ 2 

1.24e-5 

8.29e-7 

2.31e-9 

6.09e-ll 

4.78e-13 

2.18e-14 

2.02 


10" 3 

6.99e-5 

1.73e-6 

1.09e-8 

5.38e-ll 

1.17e-12 

1.59e-14 

2.22 


1 

8.04e-6 

9.05e-8 

6.80e-10 

1.39e-ll 

2.08e-13 

3.02e-15 

2.17 

(c) 

10“ 2 

3.Ole-5 

1.71e-6 

3.85e-9 

1.26e-10 

9.22e-13 

4.21e-14 

2.04 


10" 3 

1.16e-4 

4.09e-6 

2.65e-8 

6.65e-ll 

2.75e-12 

3.49e-14 

2.19 


Table 5. The L 2 errors for initial data (a)-(c) with h = 10 -5 and n(a) = X[i/ 2 ,i] (a), 
by the Laplace transform method. The notation r denotes the exponential convergence 
rate in the error ||u^ h — u(t n ) Hl^q) < Ce~ rN . 


case 

t \ N 

3 

5 

7 

9 

11 

13 

r 


1 

4.54e-6 

2.30e-7 

1.63e-9 

1.69e-ll 

2.36e-13 

8.46e-15 

2.02 

(a) 

10" 2 

6.21e-5 

1.65e-6 

3.71e-9 

1.07e-10 

7.00e-13 

2.58e-14 

2.16 


10” 3 

8.02e-4 

3.61e-6 

1 .66e-8 

4.17e-10 

3.10e-12 

6.73e-15 

2.55 


1 

4.78e-6 

4.74e-7 

2.43e-9 

3.44e-ll 

3.49e-13 

1.87e-14 

1.94 

(b) 

10” 2 

1.03e-4 

1.13e-6 

3.58e-9 

8.78e-ll 

5.04e-13 

1.93e-14 

2.24 


10~ 3 

5.12e-4 

4.79e-6 

4.95e-8 

5.23e-10 

5.15e-12 

5.58e-14 

2.29 


1 

4.79e-6 

5.61e-7 

2.75e-9 

4.07e-ll 

3.94e-13 

2.23e-14 

1.92 

(c) 

10" 2 

1.18e-4 

6.08e-7 

3.37e-9 

7.22e-ll 

2.84e-13 

8.94e-14 

2.10 


10" 3 

1.09e-4 

5.24e-6 

6 .02e-8 

5.62e-10 

5.95e-12 

1.02e-13 

2.07 


Table 6 . The L 2 errors for initial data (b) and (c) with h = 10~ 5 , //(a) = (a — 1/2) 2 
and N = 5 at small time instances t = 10 _fc , k = 4, 5, • • • ,9, by the Laplace transform 
method. 


case \ t 

10 -4 

10 -t> 

io- e 

io - 7 

10 " s 

IO" 9 

(b) 

7.05e-6 

9.39e-6 

1.58e-5 

1.75e-5 

1.81e-5 

1.82e-5 

(c) 

6.39e-6 

1.17e-5 

1.53e-5 

1.68e-5 

1.75e-5 

1.79e-5 


Table 7. The L 2 norm of the solution for initial data (a) and (c) with h = 10 -5 , 
fi(a) = (a — 1/2) 2 and N = 10 at large time instances t = 10 fc , k = 6,8, ••• ,18, 
computed by the Laplace transform method. 


case \ k 

6 

8 

10 

12 

14 

16 

18 

rate 

(a) 

3.33e-4 

2.70e-4 

2.26e-4 

1.95e-4 

1.71e-4 

1.52e-4 

1.37e-4 

1/k 

(c) 

1.06e-3 

8.54e-4 

7.17e-4 

6.17e-4 

5.41e-4 

4.82e-4 

4.34e-4 

1/k 


respectively. A first order convergence is observed for all three examples and at all time instances, showing 
the robustness of the scheme. 

To examine more closely the convergence behavior of the scheme, we consider t = 10~ fc , k = 4,5, ■ ■ • ,9,, 
and for each time instance t, divide the interval [0, t] into N = 10 subintervals. The scheme works well for 
the smooth initial data in example (a), however, it works poorly for the singular initial data in example 
(c), cf. Table 10. This behavior is predicted by Theorems 5.1 and 5.3: the error is dominated by the 
factor j for L 2 (VL) initial data. In Fig. 2, we plot the error ratio \\U^ — u(r)||/r against logr for smooth 
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10 

1/k 


Figure 1. The L 2 norm of the solution for initial data (a) and (c) at t = 10 fe , k = 
6,8, • • • , 18, by the Laplace transform method 


Table 8. The L 2 errors for initial data (a)-(c) with h = 10 4 and = (a — 1 /2) 2 , 
by the backward Euler convolution quadrature. 


case 

t \ N 

10 

20 

40 

80 

160 

320 

rate 


1 

1.82e-5 

8.78e-6 

4.31e-6 

2.12e-6 

1.Ole-6 

4.74e-7 

1.05 (1.00) 

(a) 

10" 2 

8.64e-4 

3.91e-4 

1.88e-4 

9.20e-5 

4.55e-5 

2.26e-5 

1.05 (1.00) 


10 -3 

2.17e-2 

1.10e-2 

5.51e-3 

2.76e-3 

1.38e-3 

6.92e-4 

0.99 (1.00) 


1 

4.81e-5 

2.32e-5 

1.14e-5 

5.60e-6 

2.67e-6 

1.26e-6 

1.05 (1.00) 

(b) 

10" 2 

8.11e-3 

3.87e-3 

1.88e-3 

9.29e-4 

4.61e-4 

2.30e-4 

1.03 (1.00) 


10 -3 

1.48e-2 

7.46e-3 

3.74e-3 

1.88e-3 

9.39e-4 

4.70e-4 

1.00 (1.00) 


1 

5.81e-5 

2.81e-5 

1.38e-5 

6.76e-6 

3.23e-6 

1.52e-6 

1.05 (1.00) 

(c) 

10" 2 

1.Ole-2 

4.80e-3 

2.34e-3 

1.15e-3 

5.72e-4 

2.85e-4 

1.03 (1.00) 


10 -3 

7.35e-3 

3.66e-3 

1.82e-3 

9.11e-4 

4.55e-4 

2.27e-4 

1.00 (1.00) 

Table 9. The L 2 errors 

for initial data (a)- 

(c) with h 

= 10 4 and fi(a) = 

X[i/2,i](a)> 

by the backward Euler convolution quadrature. 




case 

t \ N 

10 

20 

40 

80 

160 

320 

rate 


1 

2.20e-4 

1.06e-4 

5.20e-5 

2.58e-5 

1.28e-5 

6.40e-6 

1.02 (1.00) 

(a) 

10 -2 

1.76e-2 

8.81e-3 

4.40e-3 

2.20e-3 

1.10e-3 

5.49e-4 

1.00 (1.00) 


10 -3 

3.92e-3 

1.98e-3 

9.95e-4 

4.99e-4 

2.50e-4 

1.25e-4 

0.99 (1.00) 


1 

6.52e-4 

3.11e-4 

1.52e-4 

7.53e-5 

3.74e-5 

1.87e-5 

1.03 (1.00) 

(b) 

10" 2 

1.25e-2 

6.26e-3 

3.13e-3 

1.56e-3 

7.82e-4 

3.91e-4 

1.00 (1.00) 


10 -3 

5.76e-3 

2.88e-3 

1.44e-3 

7.18e-4 

3.59e-4 

1.79e-4 

1.00 (1.00) 


1 

7.92e-4 

3.78e-4 

1.85e-4 

9.14e-5 

4.54e-5 

2.27e-5 

1.03 (1.00) 

(c) 

10" 2 

7.40e-3 

3.71e-3 

1.86e-3 

9.28e-3 

4.64e-4 

2.32e-4 

1.00 (1.00) 


10 -3 

6.10e-3 

3.06e-3 

1.53e-3 

7.65e-4 

3.83e-4 

1.91e-4 

1.00 (1.00) 


initial data in example (a). Theorem 5.2 predicts an error estimate ||— u(t)||£ 2 (q) < crlogr^ 1 . The 
log factor l 2 (t) in Theorem 5.2 is fully confirmed by Fig. 2, and thus the corresponding error estimate is 
sharp. 
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Table 10. The L 2 errors for initial data (a) and (c) with h = 10 5 and N = 10, at time 
instances t = 10 _fe , k = 4, 5, • • • , 9, by backward Euler convolution quadrature. 


case \ t 

10 -4 

10 " a 

10 " b 

10 ' 

10 " s 

10" 9 

rate 

(a) 

2.42e-3 

1.03e-4 

7.87e-6 

7.59e-7 

7.58e-8 

7.44e-9 

1.01 (1.00) 

(c) 

7.44e-3 

5.67e-3 

4.30e-3 

3.27e-3 

2.49e-3 

1.88e-3 

0.12 (0.12) 



Figure 2. The L 2 errors for the backward Euler method for initial data (a) at small 
time instances t\ = t = 10 _fe , k = 5,6,..., 11,12. 

7. Concluding remarks 

In this work, we have presented a first rigorous numerical analysis of two fully discrete schemes (one 
based on the Laplace transform and another based on convolution quadratures) for the distributed-order 
time fractional diffusion equation with nonsmooth initial data. We have provided regularity estimates 
for the solution and developed one space semidiscrete Galerkin method and two fully discrete schemes. 
Optimal error estimates for the semidiscrete scheme were shown using an operator trick due to Fujita 
and Suzuki. The first fully discrete scheme is based on quadrature approximation of the inverse Laplace 
transform with a deformed contour of hyperbolic type, and exhibits an exponential convergence. It is 
especially suited to computing the solution at many and large time instances. The second fully discrete 
scheme is based on convolution quadrature generated by the backward Euler method, and exhibits a 
first order convergence. The sharpness of the error estimates were fully verified by extensive numerical 
experiments for both smooth and nonsmooth initial data. 

This work represents only a first step towards rigorous numerical analysis of distributed order subdif¬ 
fusion, and there are a number of avenues for further research. First, the semidiscrete and fully discrete 
schemes may be extended to the distributed order diffusion wave equation, with a nonnegative weight 
n(a) £ C[0, 2]. Second, the error estimates for the semidiscrete Galerkin scheme in the case of nonsmooth 
initial data v £ L 2 (Vl) depend on the final time T. It remains unknown how to get rid of this factor. This 
is especially important if the solution is sought for large T. Third, the assumption £ C[0,1] might 
be too restrictive and its is of much interest to relax it to \±{a) £ L°°(0,1). 

Acknowledgments 

The research of R. Lazarov and Z. Zhou have been supported in parts by NSF Grant DMS-1016525 
while that of D. Sheen by NRF-2014R1A2A1A11052429. 

References 

[1] T. M. Atanackovic, S. Pilipovic, and D. Zorica. Distributed-order fractional wave equation on a finite domain. Stress 
relaxation in a rod. Intemat. J. Engrg. Sci., 49(2):175—190, 2011. 












24 


BANGTI JIN, RAYTCHO LAZAROV, DONGWOO SHEEN, AND ZHI ZHOU 


[2] E. Bazhlekova. Completely monotone functions and some classes of fractional evolution equations, preprint, 
arXiv: 1502.04647, 2015. 

[3] E. Bazhlekova, B. Jin, R. Lazarov, and Z. Zhou. An analysis of the Rayleigh-Stokes problem for the generalized second 
grade fluid. Numer. Math., 2014. DOI-10.1007/s00211-014-0685-2 (arXiv: 1401.8049). 

[4] M. Caputo. Distributed order differential equations modelling dielectric induction and diffusion. Fract. Calc. Appl. 
Anal., 4(4):421-442, 2001. 

[5] A. V. Chechkin, R. Gorenflo, and I. M. Sokolov. Retarding subdiffusion and accelerating superdiffusion governed by 
distributed-order fractional diffusion equations. Phys. Rev. E, 66:046129, 2002. 

[6] A. V. Chechkin, R. Gorenflo, I. M. Sokolov, and V. Y. Gonchar. Distributed order time fractional diffusion equation. 
Fract. Calc. Appl. Anal., 6(3):259-279, 2003. 

[7] E. Cuesta, C. Lubich, and C. Palencia. Convolution quadrature time discretization of fractional diffusion-wave equa¬ 
tions. Math. Comp., 75(254):673—696, 2006. 

[8] K. Diethelm and N. J. Ford. Numerical analysis for distributed-order differential equations. J. Comput. Appl. Math., 
225(1):96-104, 2009. 

[9] H. Fujita and T. Suzuki. Evolution problems. In Handbook of Numerical Analysis, Vol. II, pages 789-928. North- 
Holland, Amsterdam, 1991. 

[10] I. P. Gavrilyuk, W. Hackbusch, and B. N. Khoromskij. 'H-matrix approximation for the operator exponential with 
applications. Numer. Math., 92(1):83—111, 2002. 

[11] R. Gorenflo, Y. Luchko, and M. Stojanovic. Fundamental solution of a distributed order time-fractional diffusion-wave 
equation as probability density. Fract. Calc. Appl. Anal., 16(2):297-316, 2013. 

[12] M. Hasse. The Functional Calculus for Sectorial Operators. Birkhauser-Verlag, Basel, 2006. 

[13] J. Jia, J. Peng, and K. Li. Well-posedness of abstract distributed-order fractional diffusion equations. Commun. Pure 
Appl. Anal, 13(2):605-621, 2014. 

[14] B. Jin, R. Lazarov, Y. Liu, and Z. Zhou. The Galerkin finite element method for a multi-term time-fractional diffusion 
equation. J. Comput. Phys., 281:825-843, 2015. 

[15] B. Jin, R. Lazarov, and Z. Zhou. Error estimates for a semidiscrete finite element method for fractional order parabolic 
equations. SIAM J. Numer. Anal., 51(l):445-466, 2013. 

[16] B. Jin, R. Lazarov, and Z. Zhou. On two schemes for fractional diffusion and diffusion-wave equations, preprint, 
arXiv: 1404.3800, 2014. 

[17] B. Jin, R. Lazarov, and Z. Zhou. An analysis of the LI scheme for the subdiffusion equation with nonsmooth data. 
IMA Numer. Anal., page in press, 2015. 

[18] B. Jin and W. Rundell. A tutorial on inverse problems for anomalous diffusion processes. Inverse Problems, 
31(3):035003, 40p., 2015. 

[19] J. T. Katsikadelis. Numerical solution of distributed order fractional differential equations. J. Comput. Phys., 259:11- 
22, 2014. 

[20] A. Kilbas, H. Srivastava, and J. Trujillo. Theory and Applications of Fractional Differential Equations. Elsevier, 
Amsterdam, 2006. 

[21] A. N. Kochubei. Distributed order calculus and equations of ultraslow diffusion. J. Math. Anal. Appl., 340(1):252-281, 
2008. 

[22] Z. Li, Y. Luchko, and M. Yamamoto. Asymptotic estimates of solutions to initial-boundary-value problems for dis¬ 
tributed order time-fractional diffusion equations. Frac. Calc. Appl. Anal., 17(4): 1114-1136, 2014. 

[23] Y. Lin and C. Xu. Finite difference/spectral approximations for the time-fractional diffusion equation. J. Comput. 
Phys., 225(2): 1533-1552, 2007. 

[24] M. Lopez-Fernandez, C. Palencia, and A. Schadle. A spectral order method for inverting sectorial Laplace transforms. 
SIAM J. Numer. Anal., 44(3): 1332-1350, 2006. 

[25] C. Lubich. Convolution quadrature and discretized operational calculus. I. Numer. Math., 52(2): 129-145, 1988. 

[26] C. Lubich, I. H. Sloan, and V. Thomee. Nonsmooth data error estimates for approximations of an evolution equation 
with a positive-type memory term. Math. Comp., 65(213): 1-17, 1996. 

[27] Y. Luchko. Boundary value problems for the generalized time fractional diffusion equation of distributed order. Frac. 
Cal. Appl. Anal., 12(4):409-422, 2009. 

[28] F. Mainardi, A. Mura, G. Pagnini, and R. Gorenflo. Time-fractional diffusion of distributed order. J. Vibr. Control, 
14(9-10): 1267-1290, 2008. 

[29] E. Martensen. Zur numerischen Auswertung uneigenlicher Integrale. Z. Angew. Math. Mech., 48:T83-T85, 1968. 

[30] W. McLean and K. Mustapha. Time-stepping error bounds for fractional diffusion problems with non-smooth initial 
data. J. Comput. Phys., in press. arXiv: 1405.2140, 2014. 

[31] W. McLean, I. H. Sloan, and V. Thomee. Time discretization via Laplace transformation of an integro-differential 
equation of parabolic type. Numer. Math., 102(3):497-522, 2006. 

[32] W. McLean and V. Thomee. Numerical solution via Laplace transforms of a fractional order evolution equation. J. 
Integral Equations Appl, 22(l):57-94, 2010. 

[33] M. M. Meerschaert, E. Nane, and P. Vellaisamy. Distributed-order fractional diffusions on bounded domains. J. Math. 
Anal. Appl., 379(l):216-228, 2011. 


ERROR ESTIMATES FOR DISTRIBUTED ORDER FPDE 


25 


[34] M. M. Meerschaert and H.-P. Scheffler. Stochastic model for ultraslow diffusion. Stochastic Process. Appl., 116(9): 1215- 
1235, 2006. 

[35] M. L. Morgado and M. Rebelo. Numerical approximation of distributed order reaction diffusion equations. J. Comput. 

Appl. Math., 275:216-227, 2015. 

[36] K. Mustapha, B. Abdallah, and K. M. Furati. A discontinuous Petrov-Galerkin method for time-fractional diffusion 
equations. SIAM J. Numer. Anal., 52(2):2512-2529, 2014. 

[37] I. Podlubny. Fractional Differential Equations. Academic Press, San Diego, CA, 1999. 

[38] D. Sheen, I. Sloan, and V. Thomee. A parallel method for time-discretization of parabolic problems based on contour 
integral representation and quadrature. Math. Comp., 69(229): 177-195, 2000. 

[39] D. Sheen, I. Sloan, and V. Thomee. A parallel method for time discretization of parabolic equations based on Laplace 
transformation and quadrature. IMA J. Numer. Anal., 23(2):269-299, 2003. 

[40] Y. G. Sinai. The limit behavior of a one-dimensional random walk in a random environment. Teor. Veroyatnost. i 
Primenen., 27(2):247-258, 1982. 

[41] I. M. Sokolov, A. V. Chechkin, and J. Klafter. Distributed-order fractional kinetics. Acta Phys. Polon. B, 35(4):1323- 
1341, 2004. 

[42] V. Thomee. Galerkin Finite Element Methods for Parabolic Problems, volume 25 of Springer Series in Computational 
Mathematics. Springer-Verlag, Berlin, 2006. 

[43] V. Vergara and R. Zacher. Optimal decay estimates for time-fractional and other nonlocal subdiffusion equations via 
energy methods. SIAM J. Math. Anal., 47(l):210-239, 2015. 

[44] J. A. C. Weideman and L. N. Trefethen. Parabolic and hyperbolic contours for computing the Bromwich integral. 

Math. Comp., 76(259):1341-1356, 2007. 

[45] F. Zeng, C. Li, F. Liu, and I. Turner. The use of finite difference/element approaches for solving the time-fractional 
subdiffusion equation. SIAM J. Sci. Comput., 35(6):A2976-A3000, 2013. 

[46] Y.-N. Zhang, Z.-Z. Sun, and H.-L. Liao. Finite difference methods for the time fractional diffusion equation on non- 
uniform meshes. J. Comput. Phys., 265:195-210, 2014. 

Department of Computer Science, University College London, Gower Street, London WC1E 6BT, UK 

(bangti.jin@gmail.com) 

Department of Mathematics, Texas A&M University, College Station, TX 77843-3368, USA, (lazarov@math.tamu.edu, 

ZZHOU@MATH.TAMU.EDU) 

Department of Mathematics, Seoul National University, Seoul 151-747, Korea (dongwoosheen@gmail.com) 


